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Abstract 

We consider a mean field model describing the free cooling process of a two component granular 
mixture, a generalization of so called Maxwell model. The cooling is viewed as an ordering process 
and the scaling behavior is attributed to the presence of an attractive fixed point at v = for 
the dynamics. By means of asymptotic analysis of the Boltzmann equation and of numerical 
simulations we get the following results: 

l)we establish the existence of two different partial granular temperatures, one for each compo- 
nent, which violates the Zeroth Law of Thermodynamics; 2) we obtain the scaling form of the two 
distribution functions; 3) we prove the existence of a continuous spectrum of exponents charac- 
terizing the inverse-power law decay of the tails of the velocity , which generalizes the previously 
reported value 4 for the pure model; 4) we find that the exponents depend on the composition, 
masses and restitution coefficients of the mixture; 5) we also remark that the reported distributions 
represent a dynamical realization of those predicted by the Non Extensive Statistical Mechanics, 
in spite of the fact that ours stem from a purely dynamical approach. 

PACS numbers: 02.50.Ey, 05.20.Dd, 81.05.Rm 
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INTRODUCTION 



One of the basic principles over which the Thermodynamics is built, derives from the 
simple observation that when two systems A and B are brought together so that they can 
exchange energy (or in thermal contact), after some time they reach a mutual macroscopic 
equilibrium state in the absence of energy sources, i.e. their temperatures become equal. A 
second observation is that the thermal equilibrium between a third system C and A implies 
the thermal equilibrium between B and C. These two facts, which represent the content of 
the Zeroth Law of Thermodynamics, give sense to the concept of temperature and allow us 
to build thermometers |T[ . The simplest state of aggregation of matter, namely the gaseous 
state, offers a neat example of the validity of such a principle. When in equilibrium the state 
of a mixture of different gases is characterized by a single temperature, i.e. every species 
has the same average kinetic energy per particle regardless its molecular nature. 

An interesting question related to the Zeroth Law arises quite naturally when we consider 
the behavior of Granular Gases, as assemblies of moving grains at low density are usually 
called. More in general, Granular Matter represents a relatively new and unexplored area 
of research which has attracted the attention of physicists. As ordinary matter Granular 
Matter may appear under different guises: solid-like or dense aggregates of particles where 
the motion of its constituents is negligible, liquid-like when a dense aggregate flows, and as 
a gas when the mutual distances between grains are larger that their typical linear size. Of 
course, such a classification is empirical but quite loose, since we are not allowed to employ 
concepts such as thermodynamic phases, in fact Granular Matter does not fall into the 
realm of Thermodynamics. In other words, only if we are able to show that the basic Laws 
of Thermodynamics work in the case of granular materials we can employ its results. There 
is an emerging consensus that granular materials do not achieve a proper thermodynamic 
equilibrium. Hence, even if the temperature of a granular assembly can be suitably defined, it 
does not have the properties required by standard thermodynamics, as we shall discuss below. 
We shall not digress on some interesting recent proposals aimed to define the temperature 
of granular packings 0, which is a measure of the potential energy landscape, but we 
shall confine our attention to dilute systems of granular particles, where the term granular 
temperature, defined exactly in the same way as in ordinary gases, i.e. as the average of the 
kinetic energy per particle, has been coined 0. How far pursue such an analogy?. Is the 
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granular temperature a common feature of all granular gases in mutual equilibrium, i.e. the 
quantity which has the role of determining if two system are in equilibrium with respect to 
each other with respect to exchanges of energy? If it is not so, the granular temperature 
has to be regarded just as an ensemble average of the kinetic energy of the grains, which 
may depend on the microscopic details and therefore is spoiled of one of its most useful 
characteristics. 

In the present paper we shall deal with such an issue through a simple model of granular 
mixture. To be more specific we shall describe a granular mixture by means of a model of 
the Maxwell type, i.e. a model characterized by a collision rate which is independent of the 
relative velocity of the colliding particles Q . Such a class of models has its justification in the 
case of elastic Maxwell molecules, i.e. particles interacting via a soft repulsive pair potential 
of the form r l ~ 2d . In the case of inelastic particles the constant collision rate is just a 
matter of mathematical convenience, since it does not correspond to any known microscopic 
interaction. In spite of that, it has been demonstrated by several authors fj , |§ @ , |§ @, 
that the cost of sacrificing physical realism has been widely compensated by the amount of 
exact analytical results and simplification even in the numerical treatment of the collision 
process. On the other hand the real question is: what are the physical consequences of 
the constant collision rate on our predictions? Some of these are robust with respect to 
the choice of the model others will be peculiar, as we shall discuss below. An interesting 
observation is that the homogeneous cooling states of Maxwell models and of more realistic 
models become quite similar if one measures the respective evolutions in terms of the number 
of collisions occurred up to a given instant, while they are very different in the presence of 
spatial inhomogeneities ||, [pL0| ,HTT1 , [12 . In the case of granular mixtures there is another 
aspect which is robust with respect to the choice of the model, namely the existence of a two 
temperature behavior. In fact, this has already been observed by Dufty and Garzo |13| 



in 



a study of the cooling of a mixture based on the Enskog approximation and experimentally 



by Feitosa and Menon |L4] in a mixture subject to external driving. The Maxwell model, we 
shall study, also shows such a behavior, but offers the advantage that, being considerably 
simpler than other approaches, provides a guidance to understand the global aspect of the 
problem. 

Summarizing, the choice of the model is deliberately minimalistic for the following three 
important reasons: 
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a) the simplicity of the model allows us to avoid the usual problem of the closure of the 
correlation hierarchy in statistical physics and the need of introducing approximations. In 
fact, the Boltzmann equation turns out to be exact for the model we study, so that none of 
the results we find can be ascribed to a breakdown under some peculiar conditions of some 
physically motivated approximation, such as the Boltzmann molecular chaos hypothesis; 

b) most of the work can be performed analytically and even the form of the asymptotic 
solutions of the Boltzmann equation can be estimated in closed manner with a fairly good 
accuracy; 

c) our solutions can be validated by simple simulations of the microscopic collision process. 
Perhaps it is worth to mention that as an unexpected byproduct of our study we have 

found a rich variety of behaviors of the velocity distributions as a function of the microscopic 
control parameters. The common feature of all these distributions is the presence of inverse 
power law high velocity tails. These arise quite naturally from our microscopic dynamics 
and do not need any fine tuning of the parameters, a feature somehow similar to that which 
characterize the so called Self Organized Criticality. 

Incidentally, the observed inverse-power laws for the distributions are identical to those 
predicted by the Non Extensive Statistical Mechanics (NESM) approach [|Tj| Ul6| , which 
apart from some exceptions still lacks of dynamical foundations. An important point to 
stress is the fact that the observed behavior derives from the solution of the Boltzmann 
transport equation for a system clearly inspired to the physics of granular materials and 
already employed in the past, and not the result of assumptions about the form of the 
generalized entropy, as we shall see below. 

The structure of the present paper is the following: in section II we introduce the model, 
define the notation and write the Boltzmann equation; in section III we set up the mo- 
ment expansion and characterize by a granular temperature the macroscopic state of each 
component; sections IV and V, in which we discuss how the scaling solutions for the Boltz- 
mann equation can be obtained, contain all technical details; in section VI we comment 
the connections between our results and the the NESM approach; in VII we present our 
conclusions. 



4 



DEFINITION OF THE MODEL 



We shall consider an assembly of N% particles of species 1 and N2 particles of species 2 
endowed with scalar velocities vf , with a = 1,2. The two species may have different masses, 
777.1 an d fn 2 and/or different restitution coefficients r n , r 22 , r 12 = r 21 . 

The main assumptions of the model are: 

a) the forces are pairwise and impulsive, so that their velocities change at each collision 
according to the following rule: 



«f = *-[ 1 + H^W-'*> da) 



•?=^ + i 1+ ^i^w-^ (lb) 



b) collisions involving simultaneously more than two particles are disregarded; 

c) within the spirit of mean field models all pairs are allowed to exchange momentum 
regardless of their mutual positions. 

In eqs. (|l|) the primed quantities are the post-collisional velocities and are a linear 
combination of the velocities before the collisions. The process conserves the number of 
particles of each species and the total impulse II = J2 a =i 2 Sj=i m av2- 

To study the statistical behavior of the system we consider the coupled equations for 
the velocity probability distribution functions P a (v,t), which give the probability density 
of finding a particle of species a with velocity v at the instant t. In the absence of driving 
forces we write the following two component Boltzmann equation: 

2 r 1 

d t P t (v u t) = J2 / dv 2 [-P i {v'l)P j {v'i) - Pi^Pjim)]. (2) 
j= i J a 

The system of Boltzmann equations (fj) describes the evolution of the p.d.f. of the 
velocities in the case of N\ and iV 2 going to infinity. From a rigorous point of view correlations 
can appear in a finite system even if such a mean field approach is employed, i.e. the 
Molecular Chaos assumption underlying the Boltzmann equation is no more guaranteed to 
be true for finite Ni and N 2 . The typical mechanism of failure of the Molecular Chaos is 
the re-collision process: the collisional history of two colliding particles may overlap and 
therefore the two colliding particles are somehow correlated, e.g. the two particles may have 



collided together in the recent past or may have collided with different particles that have 
collided together at previous times (ring collisions). In the simulations we have always used 
at least one million of particles: this is enough, we think, to expect a good comparison with 
analytical results. 

Eliminating the pre-collisional velocities v",v 2 ' in favor of the post-collisional velocities 
by means of eq. ([IJ) we obtain the non linear system: 



d t P x {v,t) = -P 1 (v,t)+ T ^— [duPtMPt ( 2V , (1 Tn)U ,t 



1 + r n J V 1 + n 
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(3a) 



1 + ri2 m 2 V V 1 + ri 2 



1 + r 22 J ' V 1 + r 22 

p m 1 +m 2 f mi v - - n 2 



(3b) 



+ = f 1 ' Z duP.MP mi - — ,i 

1 + r i2 mi J \ 1 + r 12 



The generalization to many components is straightforward. 

To proceed further, it is useful to employ the method of characteristic functions 
i.e. the Fourier transforms of the probability densities defined as: 



oo 

ikv 



P a (k,t) = / dve lkv P a (v,t) (4) 

J — oo 

The integral equations in Fourier space assume the more convenient form : 



(5a) 



(5b) 



dtPi(k,t) = -A(M)+pA(7iiM)A((l-7n)M) 
+ (l-p)Pi{ji 2 k,t)P 2 {l-j 12 )k,t) 

d t P 2 (k,t) = -P 2 (k,t) + (1 - P )P 2 { l22 k : t)P 2 {{l - 7 22)M) 

+ pA(72iM)A((l-72i)M) 
where p = Ni/(N\ + N 2 ), ( = m\jm 2 and 

la/3 = — - — (6a) 

7i2 = [l- I ^( 1 -7i2)] (6b) 

2 

72i = [1 - 1 + - 7i2)] (6c) 
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Interestingly, a mixture of elastic particles of unequal masses does equilibrate under the 
evolution rule (JJ), unlike the case £ = 1. In fact, the solutions of @ are Gaussians e~™^ k , 
and in real space correspond to Maxwellian velocity distributions. 

THE MOMENT EXPANSION OF THE DISTRIBUTION FUNCTIONS 

For finite inelasticity, let us tentatively assume P a (k, t) to be analytic at the origin k = 
and perform a Taylor expansion in powers of k: 

Pa(k,t) = J2 (J ^r^(t) (7) 

z — ' n\ 

n=0 

The coefficient jJ%(t) of order n of the expansion represents the n-th moment of the 
distribution P a (v, t) according to the definition eq. (|]). By substituting the Taylor expansion 
and equating coefficients of powers of k to zero one obtains a set of linear equations for the 
moments which can be systematically solved. The moments of order zero are ^(t) = 1, due 
to the normalization of P a (v, t) , while those of first order are solutions of the coupled linear 
equations: 

^i(t) = -(l-p)(l-7i 2 )[M 1) -^ 2) ] (8a) 
d t tf(t) =p(l-72i)[M 1) -/4 2) ] (8b) 

The total impulse II = N\pm\ix^ + (1 — p)m 2 /i[ 2 ' ) ] is a constant of motion and corresponds 
to the eigenvalue z\ = associated with the Galilean invariance, of the linear system (^j), 
whereas the negative eigenvalue Z2 = —(1 + r i2 ) [ mi ^ 1 +^ m2 ] > describes the exponential law 
e~ 22 ' by which two subsystems 1 and 2 in relative motion with respect to each other come 
to have the same average velocity. 

The second order cumulants represent the most important statistical indicator of the state 
of a granular gas and are usually taken as the definition of partial granular temperatures . 
In order to study their evolution one has to solve the following coupled linear equations for 
the second moments: 
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d t ^(t) = d?U i 2 ) + d { 3^ ) (9a) 
^ 2) (t) = 4 2 1 ) ^ 1) + 4 2 2 V? ) (9b) 



where the coefficients are given by: 



d&> = -1 + p[ 7 n + (1 - 7ii)1 + (1 - p) [7i 2 ]" (10a) 

4? = (l-p)[(l-7i 2 )] B (10b) 

4? = -1 + (1 -p)[ 72 n 2 + (1 - 722)1 +P[72i] n (10c) 

4i ) =P[(l-72i)] n (lOd) 

The solution of the system flj^) is a linear combination of real exponentials which can be 
expressed as: 



^(t) = A»e Alt + S tt e Aat (11) 

where Ai and A 2 are respectively the less negative and the more negative eigenvalue of the 
secular equation associated with the linear system (|^). Their computation is rather tedious 
and due to the presence of many control parameters not particularly illuminating, a part 
from simple limiting cases, for which one can extract useful information. For this reason in 
most cases we shall give their values by solving numerically the associated secular equation. 
An example of the dependence of the larger eigenvalue Ai from the control parameters is 
shown in Fig. [l|. 

For vanishing first moments the global granular temperature can be defined as: 



T g =pT 1 + (l-p)T 2 (12) 

where T a = m a [ia /2 represents the partial granular temperature of species a. 

Notice that in the limit of an elastic system (all 7^/3 = but arbitrary £ and p) the 
eigenvalue Ai = reflects the energy conservation. On the other hand, for non vanishing 
inelasticity, i.e. values of | > 70-/3 > 0, the energy is dissipated at a finite rate. Consider, for 
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instance, an arbitrary composition p, but 7 Q/ g = 7 and ( = 1, a situation which describes the 
approach to the scaling regime of two subsystems 1 and 2 initially prepared at two different 
initial granular temperatures. One finds that the "energy" eigenvalue is Ai = 27(7 — 1) 
while A2 = — (f — 7 2 ). To appreciate the role of A2 let us consider the ratio of the granular 
temperatures: 

ll = r 1 Al (13) 

Relation ( |13|) means that the resulting "thermalization" time, viz. the time spent to 
reach the homogeneous cooling regime is proportional to r = (A2 — Ai) -1 = (1 — 7)~ 2 , i.e. 
is related to the difference between the two eigenvalues and attains its maximum for the 
largest inelasticity (7 = 1/2). The asymptotic value of the temperature ratio is given by: 



T\ A (4 2 1 ) - 4 2 2 ) ) + V(4 2 1 ) - 4 2 2 ) ) 2 + 44 2 2 ) 4 2 1 ) 

T 2 C 24? { } 



It is clear from eq.(jl4j) that the two granular temperatures of the homogeneous cooling 
system are in general different (see fig. as already observed in a different model in ref. 

In the case of an equimolar mixture (p = 1/2) with equal inelasticity parameters (7^ = 
7), but arbitrary mass ratio the temperature ratio is particularly simple: 

T x _ (C 2 - 1)7 + ygj ~ C 2 )V + 4(1 - 7 )V 

T 2 ~ 2(I^K (15) 

Notice that the second moment ratio p = ^§ = ^ approaches for ( — > 00 and 
for ( — > 0. For small departures from the value ( = 1 a good approximation is: 



p~i-- — Kc-i) (16) 
1-7 

i.e. the temperature ratio increases linearly as a function of ( and deviates from the 
equipartition value 1 of an ideal gas mixture for non vanishing values of the inelasticity 
parameter 7. In Figs. ^| and [| we display the moment ratio and the temperature ratio 
respectively for different values of the mass ratio and of composition. 
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When one carries on the evaluation of higher order moments one discovers the following 
phenomenon, which is the symptom that the moment expansion is ill defined: the rescaled 
moments beyond that of order m diverge, that is to say the decay of the m-th moment is 
slower than that of the m/2-th power of the second moment. Ben-Nairn and Krapivski J7| 
found that in a pure system such a phenomenon occurs for all moments beyond the fourth. 
In other words the kurtosis of the velocity distribution diverges. This is the fingerprint of 
the presence of an inverse power-law tail in the velocity distribution function. The situation 
became clearer after the discovery of a scaling solution || as we shall discuss below. 

A simple check shows that the state of affairs remains the same in the case of a multicom- 
ponent mixture. To show that, we write the evolution equations for the fourth moments: 



dtUt = diilM + 4.2/4 + a n(/4 ) 2 + 0i2^2 4 ( 17a ) 

<9 t /4 2) = 4 4 iVi 1} + 4? /4 2) + M/if) 2 + a 2 i/i?V2 2) (17b) 

where the coefficients a a p are given by: 

an = 6p[ 7 n(l-7ii)] 2 (18a) 

a 2 2 = 6(l-p)[ 7 22(l- 722)] 2 (18b) 

a 12 = 6(l-p)[(l-7i 2 )] 2 [7i 2 )] 2 (18c) 

a 2 i = 6p[(l-72i)] 2 [72i)] 2 (18d) 



We have explored the hyper-surface ( = 1 of the five dimensional parameter space 7 a/ g, p, ( 
by randomly generating the values of the coupling constants 7 aj g and of the composition p 
and computed the eigenvalues, Ai and A 2 , with Ai > A 2 , of the secular equation for the 
fourth moments, observing that the less negative eigenvalue Aj is larger than twice the 
eigenvalue Ai. This is equivalent to say that after a transient the rescaled fourth moment 
will begin to diverge as in the one component case. Interestingly, for ( / 1 we observed 
finite rescaled fourth moments. However, this fact does not imply that for such values of the 
parameters the moment expansion of the rescaled distribution holds to all orders. It only 
means that an analogous problem appears at an higher order of the moment expansion, say 
at the m-th moment. We shall elaborate on such an issue in section V. 
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THE SCALING SOLUTION 



A clue to understand the divergence of the rescaled higher moments came recently from 
the analysis of the pure case. An exact asymptotic scaling solution of the master equation 
(0) valid in the asymptotic regime has been found || . The merit of such an explicit solution 
is to show that the characteristic function does not possess a moment expansion, because of 
the presence of a non analyticity at the origin, k = as shown by its representation: 



f(kv (t)) = e-^ k \(l + v \k\) (19) 

where lim^oo P(k, t) = f(kv (t)) and vq = Ae~ 7( ^ 7_1 ^ is the square root of the second 
moment. 

To appreciate such a result one can return to real velocity space where the distribution 
space is of the form: 



2 11 

(20) 



Mty ™ [i + (^) 2 ] 



212 



and that the moment beyond the third diverge due to the presence of the inverse-power 
v ~ 4 tail. In fact, it is a general principle that when a singularity of the function in fc-space 
approaches the real axis the behavior of the large v components of its Fourier transform 
turns from exponential to power law |17| [18 . 



An immediate check shows that in the two-component case the above formula fails to 
provide the correct solution of the master equation, with the exception of some special 
cases, which correspond to a ratio of the second moments equal to 1. 

Let us stress the importance of scaling solutions in the cooling of an inelastic gas. Their 
existence means that the distribution is fully characterized by its second moment and is 
self similar, i.e. asymptotically its shape does not change apart from a trivial rescaling. 
Moreover, it demonstrates that in the cooling problem the distribution functions do not 
look at all like the distribution functions of an elastic system, which have an infinite number 
of non diverging moments. The change induced by the presence of inelasticity is a singular 
perturbation, viz. an abrupt qualitative modification of the statistical properties of the 
systems. 
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Do scaling solutions exist in the case of binary mixtures? To answer such a question 
let us assume that for a constant cooling rate there exists a scaling solution of the form 
P a (k,t) = f a (/j,2 (t)k 2 ). This implies that the master equation for P can be converted into 
a differential equation with respect to the independent variable k with time independent 
coefficients. In the one component case it can be cast in the form: 



-akd k P(k) + P(k) - P(jk)P((l - j)k) = (21) 
with a a constant. Eq. (21) displays the presence of a regular singular point at the origin, 



k — 0. Therefore we look for a particular solution of the the form: 

y{k) = k a R{k) 

where a is the so called indicial exponent and R(k) is a function which is analytic at k — 0, 
whose value can be determined by the Frobenius method [TJJ, which provides a systematic 
tool to compute the solution as a power series of k with non integer exponents. If R(k) is 
analytic one can write: 

oo 

y{k) = k°Y, a nk n 

71=0 

However, instead of employing the general method we found more convenient to extend 
to the two component case the abridged version employed by Ernst and Brito to ascertain 
the nature of the high- velocity tails of the distributions P|,[pO[|. Their technique consists in 
replacing the small-A; Frobenius expansion with its first few terms which contain a contribu- 
tion non analytic at the origin and determining self-consistently the indicial exponent. We 
assume that the leading singularities of the Fourier transforms are of the form: 



P a (k,t) = 1 - \^\t)k 2 + {^\t)S a \k\ 2 r /2 + (^\t)) 2 Q a k 4 (22) 

with exponent a and amplitudes Q a , S a to be determined by substituting the above 
approximation into the master equation. The reason for keeping the fourth term will become 
clear below. Notice that the first two terms in eq. (|22[ ) are fixed respectively by the 
normalization condition and by the variance of the distribution. In fact, up to the second 
order eq. (p2|) is analytic and identical to the Taylor expansion (0). 
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Substituting eq..(p2[)) into eqs. (|3|) and equating the like powers of k we obtain the 
following set of coupled equations: 

Si[o t + 1 - P (r u + (i - 7ii) ct ) - (i -P)\iurM } r /2 , x 

(23a) 

-5 2 [(i-p)ii-7i 2 |](^ 2) r /2 = o 

5 2 [^+i-(i-p)( 7 f 2 +(i-7 22 r)-pi72ir](^ 2) r /2 , 1N 

(23b) 

-5 1 bii- 72 iii(4 1) r /2 = o 

To proceed further, we notice that for times much longer than the "thermalization" time 
we can safely approximate the second moments with their projection along the eigenvector 
corresponding to the larger eigenvalue. Doing so the equations Q2"3"D reduce to a homogeneous 
system for the variables S a . Only in correspondence of special values of the exponent a the 
determinant of the coefficients is zero and there exist non vanishing solutions. In other words 
we are requiring a solvability condition. The resulting indicial equation is highly non linear 
in the unknown a, but its numerical solution is straightforward. It reads: 



[_ Al + 1 - p [ 7 - + (i _ 7ll )i - (i - p) fen 

x [|a 1 + 1-(1- P )[ 72 ct 2 + (1 -722)1 

-P\l2i\ a ] -P(l -p)|l - 7i 2 |11 - 7 2 i| a = (24) 

We observe that eq. (|23j) has always two solutions a = and a = 2 for any choice of 
the parameters 7 aj/ g, C an d P- In fact, these two values correspond to the zeroth and second 
moment. The non trivial value a > 2 represents the indicial exponent associated with the 
singularity. Somehow surprisingly, we find that such a value of a is distributed continuously 
in the interval (2 < a < 3) as a function of the control parameters. In other words it means 
that the inverse-power law tails of the distribution are sensitive to the composition and to 
the nature of the interactions in the mixture. We have not found any simple dependence of 
the exponent a, on the control parameters. Nevertheless, for small asymmetry we observe 
that a seems to deviate from the exponent a = 3 of the pure system quadratically with the 
temperature ratio. 

Knowing just the first three terms of the series (^), apart from the value of the constant 
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S\ = S2 which is still at our disposal, we make the hypothesis that they represent the 
truncation of the expansion of the following characteristic function: 



2 kh 

P a (kb a ) = — ^y Kv {kb a ) (25) 

where v = a /2 and K v (z) is the modified Bessel function of the second kind of order v. 
To render the matter clearer we employ the the following Frobenius series representation 



21 



2 kb a 2 IT ^ (Ma)2n 1 (kb^yu 

T{u)^ rK ^ kha) = I» sin(irv) ^ [ T(n + 1 - v) ~ T(n + l + u) ] (26) 

considering just the first few terms 

P(kb)~l- 1 ( kba ? - ( kb «)*> r 0-- ,/ ) + ( kb Q a T(1-u) {m 

(u-l) [ 2> 1 2 } T{l + v) +[ 2 > 2T(3-u) W) 



Thus one can see that the first terms of the series (|22| ) and ( pq) have the same exponents. 
Therefore by suitably choosing the coefficients of ([22] ) we can obtain the Bessel function. 
In reality, we do so because we are led by the solution of the pure case which corresponds 
to v — 3/2. Inserting such a value and employing the following asymptotic expansion (for 

z > 0) 

K 3/2 (z) = ^f z e- z (l + z- 1 ) 

after substitution in eq. ( p5|) we obtain the solution (|i~9|). 

If we insist that the identification between the series (|22"D and the modified Bessel function 
holds for arbitrary v we obtain by Fourier transforming eq. ( p6| ) the following approximation 
to the distribution function ||22|| : 



p«ir\) = r d ^- ikvp ^) = r d ^ ikv ^-^r K ^) (28) 

v a [t) J-oo 2tt y_ 00 2tt F(v) 2 



i.e.: 



a{ v a (t) > b a ^ T{v) [l + (£)^+§ 1 ' 
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where v 2 a = 2 (l-i) represents the second moment. 

Such an approximate form clearly displays the inverse power law tails with the charac- 
teristic exponent a + 1 = 2v + 1. For a pure system since a — * 3 it reduces to the known 
solution ||. Numerically we found that such an approximation gives excellent results with 
an error which is compatible with the uncertainty of the numerical data as shown in Fig|| 
Notice that the indicial equation does not give the value of S a , but this has been fixed by 



our ansatz (25) . On the other hand we might ask how good is the ansatz. To clarify such an 
issue we substitute the expansion ( |22] ) into the scaling form of the master equation and find 
two linear inhomogeneous algebraic equations for the two unknown Qi and Q 2 . Without 
giving all the details we notice that these equations have certainly non zero solutions in 
virtue of the fact observed above, eqs. ([IT]) that the eigenvalues of the fourth moment is 
not twice the eigenvalues of the of the second moment. The numerical solution of the linear 
system gives two values of the Q's which in general are different (Qi ^ Q 2 ) for non special 
values of the parameters, hence the ansatz represents only an approximation. The conse- 
quences of this small discrepancy are probably of minor importance for the cases considered 
in the present section. The deviation of the true series from the series representation of 
the approximate solution, might manifest itself in a small asymmetry in the two rescaled 
distribution functions, and this might involve the region of small values of v. In principle 
there exist the possibility of constructing a better solution via the Frobenius method. How- 
ever, we believe that one could hardly find a closed form of the distribution functions via an 
inverse Fourier transform of simplicity comparable to that of eq (|T9|) . Finally, let us remark 
that the fairly good agreement between our approximation and the numerically computed 
solution stems from the fact that the form we propose embodies not only the three following 
basic ingredients: i) the normalization condition, ii) the correct value of the variance and iii) 
the appropriate tails of the distributions, but also some properties of the exact distribution 
which one can guess on the basis of pure physical reasoning. These properties are that 
P a {v) is symmetric w.r.t. v, is monotonically decreasing for v > 0, and is smooth. Such 
assumptions hugely restrict the class of all possible candidate distributions compatible with 
the first few terms of the expansion. In practice we match the small-/c behavior with the 
small r behavior, which is equivalent to the large /c-behavior. 
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HIGHER ORDER SINGULARITIES 



So far we have employed the approximate expansion (22) and found that on the hypersur- 
face C — 1 the indicial exponent a lies in the continuous interval 2 < a < 3 Correspondingly 
the distributions possess finite second moments, but diverging fourth moments. On the 
other hand for ( ^ 1 as anticipated in section III it is possible to observe different kinds 
of behaviors. The feature which makes the difference is what might be called the "isotopic 
effect" , i.e. induced solely by the difference of masses in a binary mixture of particles, other- 
wise identical. One observes the following phenomenon: all rescaled moments /4 /(/4 ) s ^ 2 
up to the m-th order are asymptotically finite, but diverge beyond m. 

This is indicative of solutions characterized by very large values of a. In other words, 
the allowed interval of values of a is expanded. This physically means that the exponent a 
evolves from a pronounced inverse power law to a Gaussian-like behavior as ( deviates from 
unity. Of course, the crossover from one regime to the other is determined by the inelasticity 
and by the mass ratio in the first place, but the value of the exponent does not depend on 
any simple way from the control parameters. 

Numerically we have considered different coupling parameters and obtained the results 
shown in Fig.|6| showing the trend that the exponent of the tails increases with decreasing 
inelasticity and with the difference |£ — 1|. 

In Fig. [7] we have plotted the distributions obtained with different values of (, showing 
the change of the exponent of the tails with (. We must stress the difficulty of obtaining 
clear numerical results for the exponents of the power law tails in the case of large values 
of the singularity a (or of the absolute value of the exponent itself, which is o + 1): in fact 
high inverse-power law tails need far larger statistics to be appreciated; moreover, it appears 
(numerically) that for large (negative) values of the expected power, the tail appears later in 
the v/vq domain: this means that with a finite statistics one can measure powers smaller, in 
absolute value, than those expected from the analysis of indicial equations. This phenomena 
has been observed also in ||. 

To simplify the analysis, let us consider first the master equation for large values of ( and 
equal inelasticity parameters. One sees that as ( — > oo the first equation for the distribution 
Pi is asymptotically decoupled from P 2 . 
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lim 7 12 = 1 (30a) 

£ — >oo 

lim 721 = 27 - 1 (30b) 

(—too 

(30c) 

For small inelasticity 7, the quantity 721 is negative as if the restitution coefficient for the 
collision (2 — 1) were larger than one, while 1 —712 = 0. Under such conditions, the equation 
for the distribution function of species 1 becomes effectively decoupled from species 2 as it 
can be appreciated by substituting ( |30aj )-( ]3"0~15| ) in eq. (|5[): 



dtP^k, t) = -P^k, t) + pP^-yk, t)Pi((l - 7)k, t) 

(31) 

+ {l-p)P 1 {k,t)P 2 {0,t) 

Taking into account the fact that p2(0,t) = 1 one recognizes immediately the equation 
for the p.d.f. of a pure system, i.e. of the form given by eq. (p0[). The decay rate of the 
energy is given by Ai = 27(7 — l)p. What happens to the light component? Within this 
limit the evolution equation for P 2 (k,t) looks very different from that of Pi(k,t). In fact, 
the numerical solutions indicate the existence of tails with very large values of o. How can 
we find this value of the exponent, knowing that the tails of the species 1 diverge as v ~ 4 ? 
The trick is to substitute in (5b) the expansion (l22f) and neglect the contribution stemming 



from Pi(k, t) because it corresponds to another degree of singularity. At the end one obtains 
the new indicial equation: 



27(7 - l)Py + 1 - (1 -VW 2 + (1 -lY 2 ] +P|1 - 2 7 r = (32) 

which for small values of 7 has a solution 02 ~ ^tjz^ ■ That is to say for a quasi-elastic 
system the exponent of the singularity diverges and all moments exist. To understand such 
a result we can reconsider the rule ([!]) and see that for £ — > 00 the collisions between unlike 
particles do not change neither the energy nor the velocity of the heavy species, whereas it 
changes the velocity of particles 2 in the following way: 



/(2) (2) . >/ (1) (2)\ / OQ \ 
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i.e. the species 1 acts on 2 as a stochastic noise, since t>W is randomly distributed 
according to Pi(v,t). Such a phenomenon is peculiar of the inelastic system under scrutiny 
and represents a sort of "reversed Brownian" motion in which the heavy particles act as a 
heat bath for the light particles. 



In the case of finite ( we have solved the indicial equation fl24|) and constructed the curves 
shown in Fig. §. For small values of the inelasticity 7 the surface raises steeply as ( departs 
from 1, attains a maximum and then decreases again reaching the asymptotic value a = 3 
for very large values of (, as predicted by our asymptotic analysis (see fig. |]). For larger 
values of the inelasticity such isotopic effect is less pronounced. This reentrant behavior of 
a with ( is mirrored by an analogous behavior of the rescaled moments. 

To conclude, we remark that for large values of v the coefficients of the series expan- 
sion in powers of v of order n < v are very close to the corresponding coefficients of the 
series expansion of the Gaussian, and only for n > v the power law behavior of the tails 
becomes manifest. This means that numerically it might be very difficult to detect such a 
region. Finally for large z/, which correspond to elastic systems, one recovers the Gaussian 
distribution: 

1 Y(iy + -) 1 11 -^1 

l im 1 Z,, 2 r = ^=-e ^ (34) 

The re-entrance phenomenon of the exponent a deserves a comment. In fact we have 
seen that a controls the tails of the distribution of the heavy component only. In fact, even 
for ( finite, but larger than some value ( cr one observes a bifurcation of the exponents. 
For ( > ( cr the exponent of the heavy component begin to be different from the exponent 
of the light component. Therefore, our approximation consisting in assuming the same 
indicial value for both subsystems becomes untenable. Nevertheless, we can find the correct 
values of the exponents <j\ and 02, by considering that there is no interference between the 
two singularities associated with k ai and k a ' 2 and therefore we get two decoupled indicial 
equations. The resulting scenario is quite intriguing. After the bifurcation we have two 
separate trajectories obtained by drawing the values of <j\ and o<i against (. Correspondingly 
the probability distribution functions of the two subsystems belong to two different critical 
hyper-surfaces. In the limit of 7 — > one subsystem flows to the Gaussian elastic fixed 
point, the other flows to the inelastic fixed point. 
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RELATION TO NON EXTENSIVE STATISTICAL MECHANICS 



We have discussed the cooling behavior of a simple model of granular material and found 
that the process gives rise to scaling forms of the probability distribution function. Somehow 
surprisingly, the distribution function has an inverse power law decay for large velocities of 
the form v~ a with an exponent a which depends in a non trivial fashion on the various 
control parameters and can vary in the interval 2 < a < oo in a continuous way It is an 
intriguing fact that such inverse-power law like behavior looks undeniably similar to that 
emerging in the context of the so called Nonextensive Statistical Mechanics To the 

best of our knowledge, no other model of comparable simplicity to the present yields similar 



laws [g3| as a result of a dynamical process. We stress again that none of our results has 
been derived from assumptions about the functional form of a Generalized Entropy and an 
associated maximum Entropy Principle. Our results follow from a completely independent 
approach, namely the exact treatment of inelastic collisions based on the master equation. 
A possible explanation of why the Nonextensive Statistical Mechanics and our formulas 
are so similar can be found in the non extensive properties of the energy which is a key 
assumption of NESM and a characteristic of our granular model. In addition, the mean 
field type of couplings of the present model might be a crucial ingredient to determine the 
observed inverse power law tails. 

To make an example let us take two different systems at the same temperatures and 
observe that they do not equilibrate at a temperature ratio equal to one. Repeating the 
experiment with a third system we would find a different temperature ratio for each new 
pair we can form. This reflects the non conservation of the kinetic energy; moreover the 
amount by which it is non conserved depends on the microscopic details of the two systems 
we bring together. Classically the temperature, besides being the average of the kinetic 
energy, is also the intensive thermodynamic field conjugated to the entropy and the latter in 
turn is related to the distribution function. This is not the case of the granular temperature, 
which does not possess a definition via the entropic route, but only via the kinetic route. In 
NESM, instead, one defines a generalized entropy, 



S q - 

Q 



i 

where the pi are the probabilities associated with the microstates i of the system. There- 
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fore S q does not have a unique expression as in Classical Physics, but varies with the exponent 
q, connected to a by the relation: 



2 

q = 1 + 



a + 1 

Thus S q would be a function of the various couplings of the particular system under 
scrutiny. 

Summarizing, although the NESM is capable to yield our results, it does not provide 
the value of the exponent q so that one needs an independent procedure to obtain it and 
compute the appropriate q-entropy. This non universality of the q exponent means that 
one would need a q-entropy for each particular mixture. Finally, a challenging situation for 
which we do not know the answer within the NESM, is the one inspired by the results of 
section V, where the two- components of the mixture are characterized by different power 
law tails; even within the same "experiment" one has to resort to two different q-entropies. 

This state of affairs is just the result of the absence of the Zeroth Law of Thermodynamics 
for granular systems, i.e. the fact that one needs a different thermometer for each particular 
granular mixture. 

Let us anticipate that in the case of a homogeneously heated granular system our model 
still predicts a two temperature behavior under rather general conditions, but with the 



partial temperatures not varying in time [pq| . On the other hand, the velocity distribution 
function will not show any power-law tails, but will be very similar to Gaussians and possess 
all finite moments. Thus, the power-law behavior is a peculiarity of the cooling process and 
not a necessity of non extensivity, at least as far as our model is concerned. 



CONCLUSIONS 



We have seen that the Maxwell inelastic mixture displays a continuous spectrum of expo- 
nents characterizing the inverse power-law tails of the velocity distribution functions. The 
previously found solution of the pure system p — 1 with tail exponent a + 1 = 4 can now 
be regarded as special case. In fact, a whole spectrum of exponents extending from 3 to oo 
exists. A natural question to ask is why the observed exponent a deviates from the value 
of the pure system. The naive guess is that, since in the pure case the the solution ([19]) 
shows the remarkable feature that its functional form is independent of the cooling rate, i.e. 
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of the microscopic inelasticity parameter 7 and that the tails decay invariably as v ~ 4 , the 
same should be true in the case of a mixture of particles characterized by different restitution 
coefficients. On the contrary, in the mixture case the a exponent varies with the microscopic 
parameters. Both these contrasting aspects, namely a single value of the exponent for the 
pure model and a continuous spectrum of values for the mixture recall the Renormalization 
Group theory of Critical Phenomena, where critical systems having different Hamiltonians 
( viz. different 7's in our pure case) in the neighborhood of a fixed point share the same 
exponents. That is to say that the evolution drives all systems belonging to the same critical 
(hyper) surface towards a common fixed point. This is why all pure systems exhibit identical 



asymptotic behaviors in spite of having different inelasticity parameter 7 ||24j| . On the other 
hand, mixtures correspond to systems whose trajectory lies on a hyper-surface whose points 
are attracted towards a fixed point characterized by a 7^ 3, because of the presence of other 
scaling fields. 

The model can also be viewed as a quench realized rendering suddenly inelastic an ini- 
tially elastic system. The subsequent cooling process takes the distribution from its elastic 
Gaussian fixed point to to the zero temperature fixed point. The trajectory along which this 
process occurs is characterized by peculiar values of the distributions. The scaling behavior 
is associated with the dynamics of the approach to the fixed point. 

Why the scaling solution is selected? Let us make the following observation: if we consider 
a mixture with a large number of identical components, the mixture formalism allows us 
to describe independently the evolution of different subsystems. The associated secular 
equation for M components will have M eigenvalues, M — 1 of them corresponding to the 
faster modes. The scaling solution corresponds to the slowest mode and is therefore stable, 
with respect to fluctuations. 

To conclude, the model we have investigated can be regarded as a kind of self-organized 
critical system. Due to the mean field nature of the model the fluctuations of a given 
portion of the system can influence all other elements, i.e are very long ranged. In a more 
realistic description of the inelastic interaction the scenario illustrated above will not survive 
asymptotically, because of the onset of spatial fluctuations, i.e. of correlations which tend to 
erode the high velocity tails. These tails, on the other hand, might remain observable during 
the homogeneous cooling regime but will cease after the formation of spatial gradients in 



the system [11 
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Finally, we argue that the two temperature behavior and the existence of a common 
cooling rate for the two species and of scaling solutions is a robust feature even when a more 
realistic treatment of the collision term is performed |26 . 
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FIG. 1: Maximum eigenvalue Ai of the coefficient matrix of equation (|9|) for the evolution of the 
second moment, as a function of the mass ratio Different cases are shown for various choices of 
parameters, p and m = T22 = T\i = r. 
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FIG. 2: Plot of the average kinetic energy per particle versus collision time. The two upper curves 
display the energies of system 1 and respectively. Notice the initial growth of the average energy 
of the light species 2. The lowest curve represent a system, whose temperature ratio is 1 




FIG. 3: Asymptotic ratio p = 1$ / iffi of second moments, as a function of the mass ratio (, 
for different values of the inelasticity parameter 711 = 722 = 712 = 7 = (1 — r)/2, where r is the 
restitution coefficient, p = 0.5 for all the curves. The perfectly inelastic case 7 = 1/2 corresponds 
to equal second moments for every value of the mass ratio. 
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FIG. 4: Ratio Ti(oo)/T2(oo) between asymptotic temperatures of the two components of the 
mixture, as a function of the mass ratio C, for different values of the inelasticity parameter 711 = 
722 = 712 = 7 = (1 — r )/2, where r is the restitution coefficient, p = 0.5 for all the curves. The 
perfectly elastic case is the horizontal line, which represents energy equipartition. 
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FIG. 5: Rescaled asymptotic velocity distributions from numerical simulations of the model 
with mass ratio ( = 1 and different values of the restitution coefficients. In the inset the whole 
distributions are shown, in the main graph the tails are magnified in log-log scale, and fitted with 
an inverse power law. The predictions obtained by solving the indicial equation are a = —3.86 for 
the upper curve and a = —3.21 for the lower curve. 
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FIG. 6: Power of the singularity a of the solution of the coupled master equations for the model, 
as a function of the mass ratio £ and the inelasticity parameter 711 = 722 = 712 = 7 = (— r)/2, 
with r the restitution coefficient, with p = 0.5. 
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FIG. 7: Tails of the rescaled asymptotic velocity distributions from numerical simulations of the 
model with number ratio p = 0.5 and different values of the restitution coefficients and of the mass 
ratio C- The pure case ( = 1 has the exact asymptotic solution P(x) = {2/n)(\ + x 2v J 
The Gaussian is shown guide for the eye. 
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FIG. 8: Indicial exponents a, o\ and of the singularity of the solutions of the coupled master 
equations (||) and of the singularities of the solutions of the non-coupled master equations (31) 
and (|3~2|), respectively, as functions of the mass ratio £. The squares and circles refer to the 
exponents obtained from the numerical simulations (to be compared with a\ and 02 respectively): 
the discrepancy with the theoretical predictions is due to the slow convergence, for large cr's, of 
the tails to the asymptotic value, as discussed in the text. 
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FIG. 9: Plot of the rescaled velocity distributions in the case of a large mass ratio ( = 10000. In 
the inset the two velocity distributions are shown, whereas in the figure the tails are displayed. 
The tails are characterized by a different inverse power law. 
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